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CORRELATED ELECTRON TRANSPORT IN MOLECULAR 

JUNCTIONS 

K. S. Thygesen and A. Rubio 

1 INTRODUCTION 

The dimensions of conventional silicon based electronics devices will soon be so small that 
quantum effects, such as electron tunneling and energy quantization, will begin to influence 
and eventually limit the device functionality. On the other hand molecular electronics is based 
on the idea of constructing electronic devices from bottom up using organic molecules as basic 
building blocks and in this way integrate the quantum nature of the charge carriers directly 
in the design (Cuniberti et al. 2005 and Joachim et al. 2000). Over the last decade it has 
become possible to capture individual nanostructures between metal contacts and measure 
the electrical properties of the resulting junction. The types of nanostructures vary all the 
way from a single hydrogen molecule (Smit et al. 2002) over organic molecules (Reichert et al. 
2002) to metallic chains of single atoms (Yanson et al. 1998) to carbon nanotubes suspended 
over several nanometers (Nygard et al. 2000) and inorganic nanowires and biochromophores 
(Cuniberti et al. 2005). The physics of these systems is highly non-classical showing intriguing 
phenomena such as quantized conductance, conductance oscillations, strong electron-phonon 
coupling, Kondo physics, and Coulomb blockade. For this reason a microscopic, i.e. quantum 
mechanical, understanding of nanoscale systems out of equilibrium is fundamental for the 
future development of molecular electronics. 

The theoretical description of electron transport in molecules (we often use the term 
molecule to cover a general nanostructure) represents a central challenge in computational 
nanoscience. In principle, the problem involves an open quantum system of electrons inter- 
acting with each other and the surrounding nuclei under the influence of an external bias 



voltage. Fortunately, due to the large difference in mass between electrons and nuclei, it is 
often a good approximation to regard the nuclei as classical charges fixed in their equilibrium 
positions - at least for sufficiently low temperature and bias voltage. This reduces the prob- 
lem to interacting electrons moving through the static potential created by the frozen lattice 
of nuclei. Further simplification is obtained by replacing the electron-electron interactions by 
a meanfield potential as is done in Hartree-Fock (HF) and Kohn-Sham (KS) theory. Within 
such independent-particle approximations Landauer's formula (Landauer 1970) applies, giv- 
ing the conductance as the (elastic) transmission probability for electrons at the Fermi level 
times the conductance unit Go = 2e 2 /h. Landauer's formula and, in particular, its equivalent 
formulation in terms of nonequilibrium Green's functions (NEGF) (Meir and Wingreen 1992), 
has formed the basis for almost all calculations of quantum transport in nano-scale systems. 
First-principles calculations are usually based on the KS scheme of Density Functional The- 
ory (DFT) with a local exchange- correlation functional (Taylor et al. 2001, Brandbyge et al. 
2002). These DFT transport schemes have been successfully applied to systems characterized 
by strong coupling between the molecule and the electrodes (Thygesen and Jacobsen 2005, 
Garci-Suarez et al. 2005), but systematically overestimates the conductance of weakly cou- 
pled systems (Heurich et al. 2002, Quek et al. 2007). Recently, it has been shown that the 
use of self-interaction corrected exchange-correlation functionals improves the agreement with 
experiments for such systems (Toher et al. 2005). However, such functionals contain param- 
eters which basically controls the size of the energy gap between highest occupied (HOMO) 
and lowest unoccupied (LUMO) molecular orbitals which questions the predictive power of 
such an approach. 

Apart from the problems related to the position of molecular energy levels, there are a 
number of electronic effects originating from the two-body nature of the electron-electron 
interaction, which cannot - even in principle - be described within a single-particle picture. 
These include strong correlation effects like Kondo effects and Coulomb blockade (Goldhaber 
et al. 1998, Costi et al. 1994), renormalization of molecular levels by dynamic screening 
(Neaton et al. 2006, Kubatkin et al. 2003), and life-time reduction due to quasiparticle scat- 
tering (Thygesen 2008). As we shall see in this chapter, such dynamic correlation effects can 
have a dramatic influence on the electrical properties - in particular far from equilibrium. 

In this chapter we describe how electronic correlation effects can be included in transport 
calculations using many-body perturbation theory within the Keldysh nonequilibrium Green's 
function formalism. Specifically, we use the so-called GW self-energy method (G denotes the 
Green's function and W is the screened interaction) which has been successfully applied to 
describe quasiparticle excitations in weakly correlated systems (Hybertsen and Louie 1986, 
Onida et al. 2002). To make the problem tractable, we limit the GW description to a cen- 
tral region containing the nanostructure of interest and part of the leads, while the (rest of 
the) metallic leads are treated at a mean-field level. The rationale behind this division is 
that the transport properties to a large extent are determined by the narrowest part of the 
conductor, i.e. the molecule, while the leads mainly serve as particle reservoirs (a proper 
inclusion of substrate polarization effects require that a sufficiently large part of the leads are 
included in the central region). The use of nonequilibrium many-body perturbation theory 
is only one out of several methods to include correlation effects in quantum transport. In 
another approach the density matrix is obtained from a many-body wave function and the 



non-equilibrium boundary conditions are invoked by fixing the occupation numbers of left- 
and right going states (Delaney and Greer 2004). Exact diagonalization within the molecular 
subspace has been combined with rate equations to calculate tunneling currents to first order 
in the lead- molecule coupling strength (Hettler et al. 2003). The linear response conductance 
of jellium quantum point contacts has been addressed on the basis of the Kubo formula which 
in principle allows correlation effects to be incorporated through the response function (Bokes 
et al. 2007). The time- dependent version of density functional theory has also been used as 
framework for quantum transport (Stefanucci et al. 2007, Di Ventra and Todorov 2004). This 
scheme is particularly useful for simulating transients and high frequency ac-responses and 
can in principle include correlations via non-adiabatic exchange-correlation Kernels. 

In section 2 we formulate the quantum transport problem and give a brief introduction 
to the nonequilibrium Green's function formalism. In section 3 we present the nonequilib- 
rium GW equations and discuss the important concept of conserving approximations. In 
section 4.1 we obtain an expression for the current within the NEGF formalism which holds 
for interactions in the central region. It is demonstrated, both analytically and by numerical 
examples, that a self-consistent evaluation of the GW self-energy is fundamentally important 
for nonequilibrium transport as it - in contrast to the popular non self-consistent approach - 
ensures the validity of the continuity equation. In section 5, with the aim of identifying univer- 
sal trends, we study a generic two-level model of a molecular junction. It is demonstrated how 
dynamic polarization effects renormalize the molecular levels, and a physical interpretation 
in terms of constrained total energy differences is provided. The application of a bias volt- 
age is shown to enhance the dynamic polarization effects. Moreover, quasiparticle scattering 
becomes increasingly important at larger bias leading to a significant broadening of the molec- 
ular resonances. These effects, which are all beyond the single-particle approximation, have 
large impact on the calculated IV curve. In section 6 we combine the GW-transport scheme 
with DFT (for the leads) and a Wannier function basis set, and apply it to two prototypical 
junctions, namely a benzene molecule coupled to featureless leads and a hydrogen molecule 
between infinite Pt chains, and the results are analyzed using the knowledge obtained from 
the two-level model. It is found that non self-consistent G W calculations depend crucially 
on the Go (whether it is the Hartree-Fock or Kohn-Sham Green's function). This together 
with its non conserving nature suggests that GW-transport calculations should be performed 
self-consistently. 

This chapter is a summary of recent work by the authors on incorporating many-body 
correlation effects in quantum transport, see Thygesen and Rubio 2007, Thygesen and Rubio 
2008, and Thygesen 2008. 

2 FORMALISM 

In this section we formulate the quantum transport problem and review the elements of the 
Keldysh Green's function theory needed for its solution. For more detailed introductions to 
the subject we refer to the books by Leeuwen et al. 2006 and Haug and Jauho 1998. To limit 
the technical details we specialize to the case of an orthogonal basis set and refer to Thygesen 
(2006) for a generalization to the non-orthogonal case. 



2.1 Model 



We consider a quantum conductor consisting of a central region (C) connected to left (L) and 
right (R) leads. For times t < t the three regions are decoupled from each other, each being 
in thermal equilibrium with a common temperature and chemical potentials {il,Hc, and fi R , 
respectively (see Fig. 1). At t = to the coupling between the three subsystems is switched on 
and a current starts to flow as the electrode with higher chemical potential discharges through 
the central region into the lead with lower chemical potential. Our aim is to calculate the 
steady state current which arises after the transient has died out. Notice that the duration 
of the steady state is determined by the size of the leads which we henceforth take to be infinite. 

The single-particle state space of the electrons, 7i, is spanned by the orthonormal basis set 
{(pi}. The orbitals fa are assumed to be localized such that TC can be decomposed into a sum 
of orthogonal subspaces corresponding to the division of the system into leads and central 
region, i.e. H = Hl + He + T~Lr- We will use the notation % G a to indicate that fa G H a for 
some a G {L, C, R}. The non-interacting part of the Hamiltonian of the connected system is 
written 

h = h v c L c jT (1) 

i,je cr=U 

L,C,R 

where i,j run over all basis states of the system. For a,/3 G {L,C,R}, the operator h a p 
is obtained by restricting i to region a and j to region (3 in Eq. (1). Occasionally we shall 
write h a instead of h aa . We assume that there is no direct coupling between the two leads, 
i.e. kin = km = (this condition can always be fulfilled by increasing the size of the 
central region since the basis functions are localized). We introduce a special notation for the 
"diagonal" of h, 

ho = h LL + h cc + h RR . (2) 

It is instructive to note that ho does not describe the three regions in physical isolation 
from each other, but rather the contacted system without inter-region hopping. We allow 
for interactions between electrons inside the central region. The most general form of such a 
two-body interaction is, 

V = Yl V ij,k4* c ]a> c i«' c k*- (3) 



ijkleC 



The full Hamiltonian describing the system at time t can then be written 

u(+\ f H = h + V fort<t m 
m = \H = h + V fort>t (4) 

Notice, that we use small letters for non-interacting quantities while the subscript refers to 
uncoupled quantities. At this point we shall not be concerned about the actual value of the 
matrix elements hij and V^i as this is unimportant for the general formalism discussed here. 

For times t < to each of the three subsystems are assumed to be in thermal equilibrium 
characterized by their equilibrium density matrices. For the left lead we have 

Ql = ^T exp(~P(h L - fi L N L )) (5) 



with 

Z L = Tr[exp(-P(h L - fJi L N L ))]. (6) 

Here j3 is the inverse temperature and Nl — i e i c \ a c io- is the number operator of lead 
L. qr and Zr are obtained by replacing L by R. For gc and Zq we must add V to the 
Hamiltonian in the exponential to account for the correlations. The initial state of the whole 
system is given by 

Q = QlQcQr- (7) 

If V is not included in g c we obtain the uncorrelated, or non-interacting, initial state g ni . 
Because Hq {ho) describes the contacted system in the absence of inter-region hopping, g (g n i) 
do not describe the three regions in physical isolation. In other words the three regions are 
only decoupled at the dynamic level for times t < t . 




Left lead 
(L) 



Central region 
(C) 



Right lead 

(R) 



A 



CD 
^_ 

CD 
LU 




4^ 



4^ 
4^ 



Mr 



Figure 1: Before the coupling between the leads and central region is established, the 
three subsystems are in equilibrium with chemical potentials fix, and /iq, respectively. 
Reprinted with permission from Thygesen and Rubio, Phys. Rev. B 77, 115333 (2008). 
Copyright 2008 by the American Physical Society. 



2.2 Contour-ordered Green's function 

To treat nonequilibrium problems it is useful to extend the time-propagation operator from 
the real time axis to the so-called Keldysh contour, Cj, depicted in Fig. 2(a). The contour 
is defined in the complex plane and runs along the real axis from to to infinity, then back to 
t and vertically down to to — ifl. When r and r' denote points on the Keldysh contour, the 
generalized time-propagation operator is defined by 



U(t',t) =Te- i $T Affl{f \ 



(8) 



where T orders a product of operators according to their time argument on the contour 
(later times further to the left). The integral is taken along Cj. So far we have defined the 
Hamiltonian H{r) only for r on the real axis. For Eq. (8) to make sense we must extend H(t) 
to the vertical branch of Ci which we will do in a moment. The contour-ordered single-particle 
GF relevant for our transport problem is defined by 

GirjAr, r') = -iTrigTicH^c^Ar')}}, (9) 

where e.g. CH,i a ( T ) — U(to,T)c ia U(T,t ) and g is the state of the system at time t . Notice 
that when evaluating cn,ia{ T ) f° r T on the real axis it does not matter whether r is chosen 
on the upper or lower part of the contour. Notice also that when r and r' are both in the 
vertical branch, G is nothing but the Matsubara GF known from the equilibrium theory. For 
later reference we also note that a non-equilibrium GF is completely defined once (i) the 
Hamiltonian governing the dynamics and (ii) the initial state have been specified. Since the 
Hamiltonian contains no spin-flip processes the GF is diagonal in spin space, i.e. G ia j a r = 
Saa'Gij and we therefore suppress the spin-dependence in the following. 
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Figure 2: (a) The Keldysh contour C\. The dynamics of the system is governed by the 
Hamiltonian H(t) on the horizontal branch, while the initial state is defined by H(t) on the 
vertical branch, (b) The Keldysh contour C, used when correlations in the initial state are 
neglected. 



If we define H(t) along the vertical part of the contour to be 



H{r)= (ha-VaN a ) + V, (10) 

a=L,C,R 

we see that U(t — i(3,t ) = Zg. We use this observation to write 

Gij(T,T) = -l ~ ' . (11) 

Tr{Te Jc ' ( '} 

Here the time- argument of q(t) and cj(r') only serves to identify their position in the contour- 
ordering, i.e. they are not in the Heisenberg picture. In order to obtain an expansion of 
Gijij, t') in powers of V, we switch to the interaction picture defined by the non-interacting 
Hamiltonian h{r) = H{r) — V . In this picture we have 

r f n Tr{g m T[e-^^(.) Cfe;t(r)c t , (r0]} 

Gij(r, t = -% , 12 

Tr{^Te- l ^ d ^ W } 



where the time-dependence of the operators is governed by the evolution operator in Eq. (8) 
with H replaced by h. The density matrix g ni is given by 

2 _ exp(-/3 J2 a (h a ~ H a N a )) ^ 



which differs from g in that it does not contain interactions in the central region. From the 
identity 

Zg = Z o g m Te~ l ^r 0d ^ (14) 

it is clear that the integration along the vertical branch of Ci in Eq. (12) accounts for the 
correlations in the initial state of region C. While it must be expected that the presence 
of initial correlations will influence the transient behavior of the current, it seems plausible 
that they will be washed out over time such that the steady state current will not depend on 
whether or not correlations are present in the initial state. In practice the neglect of initial 
correlations is a major simplification which allows us to work entirely on the real axis avoiding 
any reference to the imaginary time. For these reasons we shall neglect initial correlations in 
the rest of this paper. The GF can then be written 

Gj, (r, t ) = -i 1 1 , (15) 

Tr{£ ra Te-/e<W) } 

where the contour C is depicted in Fig. 2(b). Eq. (15) is the starting point for a systematic 
series expansion of Gij in powers of V. Since g ni is a mixed state of Slater determinants and 
the time-evolution is given by the non-interacting Hamiltonian, h, Wick's theorem applies 
and leads to the standard Feynman rules with the exception that all time integrals are along 
the contour C and all Green's functions are contour-ordered. The Feynman diagrams should 
be constructed using the Green's function defined by g n % and h, 

g tJ (r, t') = -iTr{^T[c M (r)4 ; .(r')]} (16) 

which describes the non-interacting electrons in the contacted system. 



The diagrammatic expansion leads to the identification of a self-energy, E, as the sum of 
all irreducible diagrams with no external vertices. The GF is related to the self-energy and 
the non-interacting GF through Dyson's equation 

G(r, r') = g(r, r') + J dr 1 dr 2 y(r, n)E(n, r 2 )G(r 2 , r'), (17) 

where matrix multiplication is implied. As we will see in Sec. 4.2, only the Green's function 
of the central region is needed for the calculation of the current, and we can therefore focus 
on the central-region submatrix of G. Since the interactions are limited to the central region, 
the self-energy matrix, Ejj, will be non-zero only when both i,j G C, and it should therefore 
be safe to use the notation E instead of E c . Restricting Eq. (17) to the central region we 
have 

G c (r,r') = g c (T,T>) 

+ Jdr^gcfanMn^Gcfoy). (18) 



The free propagator, gc{r,T f ), is a non-equilibrium GF because g n i is not a stationary state 
of h, i.e. [g ni , h] ^ 0. It is, however, not difficult to show that gc satisfies the following Dyson 
equation 

gc(r,r') = #o,c(t,t') + J dT 1 dT 2 go,c{ T i T i) 

[E L (n,r 2 ) + S iJ (r 1 ,r 2 )]( 7c (r 2 ,r / ), (19) 

where g is the equilibrium GF defined by (5 n j and /i . The coupling self-energy due to lead 
a = L, R is given by 

S a (r, t') = h Ca go,a( T i T ')h a c- (20) 

Notice the slight abuse of notation: S a is not the aa submatrix of S. In fact Si, and E# are 
both matrices in the central region indices only. Combining Eqs. (18) and (19) we can write 
the Dyson equation for Gc, 

G c (t,t') = 9o,c(t,t') 

+ J dT 1 dT 2 go ! c(r,T 1 )T, tot (T 1 ,T 2 )Gc(T2,T'), (21) 

in terms of the equilibrium propagator of the non- interacting, uncoupled system, g , and the 
total self-energy 

Z tot = S + S L + (22) 

The total self-energy describes the combined effect of the interactions and the coupling to the 
leads. 



2.3 Real-time Green's functions 



In order to evaluate expectation values of single-particle observables we need the real-time 
correlation functions. We work with two correlation functions, also called the lesser and 
greater GFs, defined as 

Gg(M') = iTr{QnC* Hij (!f)c Hti (t)} (23) 
G>.(M') = -^{QniCHA^H,^')}. (24) 

Again, the use of Q n i instead of g amounts to neglecting initial correlations. Two other 
important real-time GFs are the retarded and advanced GFs defined by 

G^.(M') = eit-t'^M-Gf^f)) (25) 
G?.(f,f) = 9{t'-t){Gf 3 {tX)-G>^t')). (26) 

The four GFs are related through 

G> - G< = G r - G a . (27) 

The lesser and greater GFs are just special cases of the contour-ordered GF. For example 
G < (t, t') = G(t, t') when r = t is on the upper branch of C and r' = t' is on the lower branch. 
This can be used to derive a set of rules, sometimes referred to as the Langreth rules, for con- 
verting expressions involving contour-ordered quantities into equivalent expressions involving 



real-time quantities. We shall not list the conversion rules here, but refer to Haug and Jauho 
1998 (no initial correlations) and Leeuwen et al. 2006 (including initial correlations). The 
usual procedure in non-equilibrium is then to derive the relevant equations on the contour 
using the standard diagrammatic techniques, and subsequently converting these equations to 
real time by means of the Langreth rules. An example of this procedure is given in Sec. 3.2 
where the non-equilibrium GW equations are derived. 

2.3.1 Equilibrium 

In equilibrium, the real-time GFs depend only on the time difference t' — t. Fourier transform- 
ing with respect to this time difference then brings out the spectral properties of the system. 
In particular the spectral function 

A(u) = i[G r {u) - G a {u)] = i[G>(c<;) - G < (u)\ (28) 

shows peaks at the quasiparticle energies of the system, i.e. at the energies E n (N +1) — E m (N) 
and E m {N) — E n (N — 1), where n,m denote energy levels and iV the number of electrons. 
We thus see that the single-particle GF carries information about the electron addition- and 
removal energies. Clearly these are the types of excitations which are relevant in a transport 
situation where electrons are continuously added to and removed from the central region. In 
section 5.2 we discuss how many-body polarization effects renormalize the spectral function 
of a molecule adsorbed at a metal surface. In equilibrium we furthermore have the important 
fluctuation-dissipation theorem which relates the correlation functions to the spectral function 
and the Fermi-Dirac distribution function, /, 

G < (u) = if{uj)A{uj) (29) 
G > {u) = -i(l-f(u))A(u). (30) 

The fluctuation-dissipation theorem follows from the Lehman representation which no longer 
holds out of equilibrium, and as a consequence one has to work explicitly with the correlation 
functions in non-equilibrium situations. 

2.3.2 Non-equilibrium steady state 

We shall make the plausible assumption that in steady state, all the real-time GFs depend 
only on the time-difference t' — t. Moreover, we take the limit t — > — oo. This will allow us to 
use the Fourier transform to turn convolutions in real time into products in frequency space. 
Applying the Langreth conversion rules to the Dyson equation (21), and Fourier transforming 
with respect to t' — t then leads to the following expression for the retarded GF of the central 
region 

CT c (u) = gl c (u) + gld^MG^)- (31) 
This equation can be inverted to yield the closed form 

CT c {u) = [(" + We -h c - E£M - TT R {u) - H'M]" 1 . (32) 
The equation for G a is obtained by replacing r by a and r\ by — i] or, alternatively, from 



= [G r c (oj)]\ 



(33) 



which follows from the assumption that the GFs depend on time differences only. For the 
correlation functions the conversion rules leads to the expression 

= G^E^HQM + A</» (34) 

where 

A</» = [I c + G r c (u) Y7 tot (u)]g^ / J > (u) [I c + ^ ot (co)G a c (co)}. (35) 

Using that E^" = (g^c) 1 ~ (Cp") -1 together with the equilibrium relations g^ c = —f(oo — 
Vc)[g r ,c ~ 9o,c\ and 9o,c = -(/(^ - Vc) - l)[^S,c ~ 9o,c\> we find that 

A<H = 2t V f(uj-nc)G r c (uj)G a c (u) (36) 
A» = 2 lf| [/( W -^)-l]Q( W )QH. (37) 

If the product G r G a is independent of rj we can conclude that A — > in the relevant limit 
of small r). It can in fact be shown that A always vanishes for interacting electrons while 
for non-interacting electrons A vanishes only when there are no bound states (Thygesen and 
Rubio 2008). 

2.3.3 Non-interacting electrons and lead self-energy 

In the special case of non-interacting electrons, the retarded and advanced GFs are indepen- 
dent of the initial state of the system, i.e. of the g entering the GF. Moreover, if the dynamics 
is governed by a time-independent Hamiltonian, g r and g a depend only on the time difference 
t' — t (even if the initial state is not a stationary state). In this case the Fourier transform of 
the retarded and advanced GFs with respect to t' — t equals the resolvent of the Hamiltonian 
matrix h, 

g r ' a (u) = [(u±ir t )I-h]- 1 , (38) 

where I is the identity matrix and i] is a positive infinitesimal. In our transport problem, the 
block-diagonal structure of h allows us to obtain the non-interacting GF of the uncoupled 
system by inverting each block separately, 

g r o / *(u) = [(u±iv)I-h a ]- 1 (39) 

for a G {L,C,R}. Now it is in fact easy to show that the central region component of g'^ a 

ScltlsflGS 

gt{u) = [(co ± ir,)I -h c - (uo) - ^»]-\ (40) 
with the retarded/ advanced coupling self-energies given by 

Z^ico) = h Ca gZ{uo)h aC . (41) 

Eqs. (40) and (41) give the retarded and advanced components of Eqs. (19) and (20), re- 
spectively. Notice that Sj/° depends on the applied bias voltage through g ,» because the 
self-consistent field in the leads follow the chemical potential to ensure charge neutrality as 
sketched in the lower panel of Fig. 1. Assuming a symmetrically applied bias (/j>l/r = ep±V/2) 
we have 

(V; W ) = E2 /a (0;lo-V/2), (42) 
with a similar relation for S R . Since g is an equilibrium GF its lesser and greater compo- 
nents, and thus also Sq^, follow from the fluctuation-dissipation theorem. In contrast, the 
correlation functions derived from g, which is a non-equilibrium function, must be calculated 
using the Keldysh equation (34). 



3 MANY-BODY SELF-ENERGY 



In this section we discuss two specific approximations to the many-body self-energy E intro- 
duced in Eq. (17), namely the GW and second Born (2B) approximations. Strictly speaking 
the E of Eq. (17) contains the full effect of the interactions whereas the GW and 2B self- 
energies only describe exchange and correlation, i.e. they do not include the Hartree potential. 
The GW self-energy is obtained by summing an infinite set of Feynman diagrams - one dia- 
gram at each order of the interaction - while the 2B approximation is exact to second order 
in the interaction (if performed self-consistently). In section 3.1 we introduce an effective 
interaction which leads to a particularly simple form of the self-energy equations and at the 
same time provides a means for reducing self-interaction errors in higher order diagrammatic 
expansions. In sections 3.2 and 3.3 we discuss the nonequilibrium GW and 2B equations 
using the effective interaction. 

3.1 Effective interaction 

The direct use of the full interaction Eq. (3) in a diagrammatic expansion is problematic as 
it introduces frequency dependent, four-index quantities, which quickly becomes difficult to 
store and handle numerically. For this reason we consider instead the effective interaction 
defined by 

VcS = Vi<T,3<7' C lv C la' C 3<T' C i<ri ( 43 ) 

ij,aa' 

where 

This expression follows by restricting the sum in the full interaction Eq. (3) to terms of the 
form Vij^ijCiaCjaiCjo'Cio and VijjiC icr Cj a Cj a Ci a . 

The effective interaction is local in orbital space, i.e. it is a two-point function instead 
of a four-point function and thus resembles the real-space representation. Note, however, 
that in contrast to the real-space representation Vi a j a / is spin-dependent. In particular the 
self-interactions, V ia;ia , are zero by construction and consequently self-interaction (in the or- 
bital basis) is avoided to all orders in a perturbation expansion in powers of V. Since the 
off-diagonal elements (i ^ j) of the exchange integrals are small for localized basis func- 
tions, the main effect of the second term in Eq. (44) is to cancel the self-interaction in the 
first term. 

It is not straightforward to anticipate the quality of a many-body calculation based on 
the effective interaction (43) as compared to the full interaction (3). Clearly, if we include all 
Feynman diagrams in E, we obtain the exact result when the full interaction (3) is used, while 
the use of the effective interaction (43) would yield an approximate result. The quality of this 
approximate result would then depend on the basis set, becoming better the more localized 
the basis functions and equal to the exact result in the limit of completely localized delta 
functions where only the direct Coulomb integrals will be non-zero. However, when only 
a subset of all diagrams are included in E the situation is different: In the GW approximation 
for instance, only one diagram per order (in V") is included, and thus cancellation of self- 
interaction does not occur when the full interaction is used. On the other hand the effective 



interaction (44) is self-interaction free (in the orbital basis) by construction. The situation can 
be understood by considering the lowest order case. There are only two first order diagrams 
- the Hartree and exchange diagrams - and each cancel the self-interaction in the other. 
More generally, the presence of self-interaction in an incomplete perturbation expansion can 
be seen as a violation of identities of the form {-\c\ a , ■ ■ ■ Ci a Ci a ■ • • Cj a />\-) = 0, when not all 
Wick contractions are evaluated. Such expectation values will correctly vanish when the 
effective interaction is used because the prefactor of the Ci a Ci a operator, \4t,mt, is zero. The 
effect of self-interaction errors in (non-self consistent) GW calculations was recently studied 
for a hydrogen atom (Nelson et al. 2007). For a quantitative estimate of the quality of GW 
calculations based on V e s we refer to Thygesen and Rubio 2008. 



3.2 Non-equilibrium GW self-energy 



It is useful to split the full many-body self-energy into its Hartree and exchange-correlation 
parts 

E(r,r') = S fc (r,r , ) + S TC (r,r'). (45) 

The Hartree term is local in time and can be written £/i(r, r') = H h (r)Sc(T, r') where 5q is a 
delta function on the Keldysh contour. The xc- self-energy is nonlocal in time and contains all 
the complicated many-body effects. In the GW approximation the xc- self-energy is written 
as a product of the Green's function, G, and the screened interaction, W. Usually W is 
calculated in the random-phase approximation (RPA), and it has been found that improving 
W beyond RPA has little effect on the resulting GW self-energy (Verdozzi et al. 1995). 
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Figure 3: Feynman diagrams for the GW self-energy in terms of the Green's function (G), 
the screened interaction (W), and the polarization bubble (P). All quantities are functions 
of two Keldysh times, and two basis function indices. Integration/summation over internal 
times/indices is implied. 



With the effective interaction (43) the screened interaction, W, and the polarization func- 
tion, P, are reduced from four- to two-index functions. For notational simplicity we absorb 
the spin index into the orbital index, i.e. (icr) — > i (but we do not neglect it). The expression 



for the contour-ordered GW self-energy in terms of contour-ordered quantities then read 



S G ^(r,r') = iG 13 (t,t' + )W 13 {t,t') (46) 

W l3 iT,T>) = %6 C (r,T') + V f dT 1 V ik P kl (T,T 1 )W lj (T 1 ,T') (47) 

kl Jc 

P 3 (t,t') = -iG ij (T,T , )G ji (r',T). (48) 



It is important to notice that in contrast to the conventional real-space representation of the 
GW self-energy, the spin-dependence cannot be neglected when the effective interaction is 
used. The reason for this is that V is spin-dependent and consequently the spin off-diagonal 
elements of W will influence the spin-diagonal elements of G, £, and P. A diagrammatic rep- 
resentation of the GW approximation is shown in Fig. 4. As they stand, equations (46)- (48) 
involve quantities of the whole system (leads and central region). However, since Vij is non- 
zero only when i,j G C, it follows from Eq. (47), that W and hence £ also have this structure. 
Consequently, the subscript C can be directly attached to each quantity in Eqs. (46)-(48), 
however, for the sake of generality and notational simplicity we shall not do so at this point. 
It is, however, important to realize that the GF appearing in the GW equations includes the 
self-energy due to the leads. 

Using the Langreth conversion rules the retarded and lesser GW self-energies become (on 
the time axis), 

= *Gl ] (t)W>(t)+ i G^(t)Wl j (t) (49) 
ESfe,-(f) = iGf^m^it), (50) 

where we have used the variable t instead of the time difference t' — t. For the screened 
interaction we obtain (in frequency space), 

W r (u) = V[I - P'i^V}- 1 (51) 
W <!> {u) = W'^P^^W^u;). (52) 

where all quantities are matrices in the indices i, a and matrix multiplication is implied. 
Notice that the spin off-diagonal part of V will affect the spin-diagonal part of W r through 
the matrix inversion. Finally, the real-time components of the irreducible polarization become 

P[ 3 (t) = - iG \ 3 it)G<{-t) - iG< 3 \t)Gl{-t) (53) 
P5 !> it) = -iG^WG^i-t). (54) 

From their definitions it is clear that both the polarization and the screened interaction obey 
the relations P^(ou) = PJ^—u) and W£j(u) = W^—uj), while for the self-energy and GFs we 
have Tiq W {u) = Y7q W (uo)^ and G a (u) = G r {uy . In addition all quantities fullfill the general 
identity X > — X < = X r — X a . We mention that the nonequilibrium GW approximation 
has previously been used to study band-gap renormalization in excited GaAs (Spataru et al. 
2004). 

In deriving Eqs. (51,52) we have made use of the conversion rules 8^ > {t,t 1 ') = and 
5^ a (t,t') = 5(t — t'). With these definitions the applicability of the Langreth rules can be 



extended to functions containing delta functions on the contour. Notice, however, that with 
these definitions relation (25) does not hold for the delta function. The reason why the delta 
function requires a separate treatment is that the standard Langreth rules are derived under 
the assumption that all functions on the contour are well behaved, e.g. not containing delta 
functions. We stress that no spin symmetry has been assumed in the above GW equations. 
Indeed by reintroducing the spin index, i.e. % — > (io~) and j — > (jcr'), it is clear that spin- 
polarized calculations can be performed by treating G^ and G±± independently. 

Within the GW approximation the full interaction self-energy is given by 

S(r, t') = E h (r, r') + E GW ,(r, r'), (55) 

where the GW self-energy can be further split into an static exchange and a dynamic corre- 
lation part, 

^cw(r, t') = T, x 5 c {r, t') + S corr (r, t'). (56) 
Due to the static nature of and S x we have 

= S</> = 0. (57) 

The retarded components of the Hartree and exchange self-energies become constant in fre- 
quency space, and we have (note that for S fc and £ x we do not use the effective interaction 
(43)) 

^ = ~ i ^2 G ki( t = ®) v ik,ji (58) 

kl 

ki 

Due to (57), it is clear that Eq. (50) yields the lesser /greater components of S corr . Since 
^cott( t , t ') does not contain delta functions its retarded component can be obtained from the 
relation, 

KorM = e(-t)[j:> w (t)-j:< w (t)]. (60) 

The separate calculation of Y7 X and S^orr fr° m Eqs. (59), (60) as opposed to calculating their 
sum directly from Eq. (49), has two advantages: (i) It allows us to treat T, x , which is the dom- 
inant contribution to T^cw, at a higher level of accuracy than S corr . (ii) We avoid numerical 
operations involving G r and W r in the time domain. For more detailed discussions of these 
points see appendices A and E of Thygesen and Rubio 2008. 

3.3 Non-equilibrium second Born approximation 

When screening and/or strong correlation effects are less important, as e.g. in the case of 
small, isolated molecules, the higher-order terms of the GW approximation are small and it 
is more important to include all second order diagrams (Stan et al. 2006). The full second 
order approximation, often referred to as the second Born approximation (2B), is shown dia- 
grammatically in Fig. 4. As we will use the 2B for comparison with the GW results we state 
the relevant expressions here for completeness. 



On the contour the 2B self-energy reads (with the effective interaction (43)) 



kl 

- G ^ r')G kl (r', r)G w (r, r')V a V ]k 

kl 

(61) 

Notice that the first term in S 2 b is simply the second order term of the GW self-energy. From 
Eq. (61) it is easy to obtain the lesser /greater self-energies, 

kl 

- Y. G * > ( t ) G T{-wfJ > miv jk , 

kl 

where t has been used instead of the time difference t — t'. Since these second order con- 
tributions do not contain delta functions of the time variable, we can obtain the retarded 
self-energy directly from the Kramers-Kronig relation 

Z r 2B (t) = 6(-t)[Z> B (t)-Z< B (t)]. (62) 

4 CURRENT 

As mentioned earlier, the current flowing between the electrodes can be calculated from the 
Green's function of the central region. After a short introduction to the concept of conserving 
approximations in section 4.1, we present the relevant formulas for evaluating the current 
under non-equilibrium conditions. In section 4.3 we then derive a condition on the interaction 
self-energy which ensures charge conservation in the sense that the current is the same at the 
left and right interface between the central region and the leads. Finally, in section 4.4 we show 
that this condition is always fulfilled for the so-called conserving, or ^-derivable, self-energies. 

4.1 Conserving approximations 

Once a self-energy has been obtained the GF follows from the Dyson and Keldysh equations 
(32), (34). Any single-particle observable, such as the current or the density, can then be 
calculated from the GF. An important question is then whether the calculated quantities 
obey the fundamental conservation laws like charge- and energy conservation. In the context 
of modeling electron transport, the condition for charge conservation as expressed by the 
continuity equation 

^n(r,t) = -V-j(r,t) (63) 

is obviously of particular interest. 

Baym and Kadanoff (Baym and Kadanoff 1961, Baym 1962) showed that there exists a 
deep connection between the self-energy and the validity of the conservation laws. Precisely, 
any self-energy which can be written as a functional derivative, = 5&[G]/5G, where 



belongs to a certain class of functionals of G, produces a GF which automatically fulfills 
the basic conservation laws. A self-energy which can be obtained in this way is called <E>- 
derivable. Since a ^-derivable self-energy depends on G, the Dyson equation must be solved 
self-consistently. The exact can be obtained by summing over all skeleton diagrams 

constructed using the full G as propagator. By a skeleton diagram we mean a closed diagram, 
i.e. no external vertices, containing no self-energy insertions. Practical approximations are 
then obtained by including only a subset of skeleton diagrams. Two examples of such approx- 
imations are provided by the GW and second Born $-functionals and associated self-energies 
which are illustrated in Fig. 4. Another example is provided by the Hartree-Fock approxima- 
tion. Solving the Dyson equation self-consistently with one of these self-energies thus defines 
a conserving approximation in the sense of Baym. 




Figure 4: The GW and second Born self-energies, Tiqw and S2B, can be obtained as functional 
derivatives of their respective <£>-functionals, ^gw^] and $2b[G]. Straight lines represent the 
full Green's function, G, i.e. the Green's function in the presence of coupling to the leads 
and interactions. Wiggly lines represent the interactions. Reprinted with permission from 
Thygesen and Rubio, Phys. Rev. B 77, 115333 (2008). Copyright 2008 by the American 
Physical Society. 

The validity of the conservation laws for ^-derivable self-energies follows from the invari- 
ance of $ under certain transformations of the Green's function. For example it follows from 
the closed diagramatic structure of <3> that the transformation 

G(vt, rV) -> e lA(rT) G(rr, rY)e- iA(rV) , (64) 



where A is any scalar function, leaves $[G] unchanged. Using the compact notation (r 1; r±) = 
1, the change in <£> when the GF is changed by SG can be written as 5$ = J dld2E(l, 2)5G(2, 1 + ) 
0, where we have used that E = 5$[G]/5G. To first order in A we then have 

5$ = iy"dld2E(l,2)[A(2) - A(1)]G(2,1 + ) 

= i y'dld2[S(l,2)G'(2,l + ) -G(1,2 + )S(2,1)]A(1). 

Since this hold for all A (by a scaling argument) we conclude that 

d2[E(l, 2)G(2, 1+) - G(l, 2 + )E(2, 1)] = 0. (65) 



This condition ensures the validity of the continuity equation (on the contour) at any point 
in space. In the following sections we derive and discuss this result in the framework of the 
transport model introduced in Section 2. 



4.2 Current formula 



When the coupling between the leads and the central region is established, a current will 
start to flow. The particle current from lead a into the central region is given by the time 
derivative of the number operator of lead a (Meir and Wingreen 1992), 

I a (t) = -i([H,N a ](t)) 

= iJ2G<(t,t)hin-h ni G<(t,t). (66) 

n£C 

A simple diagrammatic argument shows that (i G a, n G C) 

G ni (T,T') = ^2 & T l G nm{.T,T l )h mj g 0>ji (T 1 ,T') 

mec 

G in (T,r') = ^2 dT i9o,ij( r , T i)h jm G mn ( y T 1 ,T l ). 

meC 

Using Eq (20) we notice that Eq. (66) can be written as iTr[A < (t, t)] when A is defined as in 
Eq. (92) with B = G c and C = E a . From the general result (93) it then follows that 

Ia = I ^K^^M " £>(c)G<M], (67) 

where matrix multiplication is implied. By writing / = (I L — Ir)/2 we obtain an expression 
which is symmetric in the L, R indices, 

I=^f Tr[(r L - Y R )G< + (f L T L - f R T R )(G r c - G a c )]duj (68) 

where we have suppressed the u> dependence and introduced the coupling strength of lead a, 
T a = i[£ r a — E°]. We notice that when interactions are present, the integrals in Eqs. (68) and 
(67) will have contributions outside the bias window, /il < u < fi R , because the conduction 
electrons can gain or lose energy by interacting with other electrons in the central region. 



4.3 Charge conservation 

Due to charge conservation in the steady state we expect that II = —Ir = I, i.e. the current 
flowing from the left lead to the molecule is the negative of the current flowing from the right 
lead to the molecule. Below we derive a condition for this specific form of particle conservation. 

From Eq. (67) the difference between the currents at the left and right interface, AI = 
Il + Ir, is given by 

AI = I ^ + S ^ " ( ^ + ^ >r) °^ (69) 

To obtain a condition for AI = in terms of E we start by proving the general identity 

J ^Tr[S<,HG>H - S>,HG<H] = 0. (70) 

To prove this, we insert G < l > = G^T^fJ^G^, + A < / > (from Eq. (34)) in the left hand side of 
Eq. (70). This results in two terms involving G r Ti^ > G a and two terms involving A < / > . The 
first two terms contribute by 

/ ^ Tr \-^tG r ^ ot G a - Xf ot G%< t G«] . (71) 

Inserting E> oi = E< t + (G°)- 1 -(G r )" 1 ( we use that (G' a ) _1 -(G' r ) _1 = S^-E^ = £> t -£< t ) 
in this expression and using the cyclic invariance of the trace, it is straightforward to show 
that Eq. (71) vanishes. The two terms involving A < / > contribute to the left hand side of 
Eq. (70) by 

J ^Tr[E< t (u;)A>) - £>,HA<H]. (72) 

As mentioned in section 2.3.2, A < and A > are always zero when interactions are present. In 
the case of non-interacting electrons we have S^,{ > = '^l > + ^ < r > , which vanish outside the 
band width the leads. On the other hand A < / > is only non-zero at energies corresponding to 
bound states, i.e. states lying outside the bands, and thus we conclude that the term (72) is 
always zero. 

From Eqs. (69) and (70) it then follows that 

A/ = J ^Tr[S<HG>H - E>( W )G<H]. (73) 

We notice that without any interactions particle conservation in the sense AI = is trivially 
fulfilled since S = 0. When interactions are present, particle conservation depends on the 
specific approximation used for the interaction self-energy, S. 

4.4 Charge conservation from ^-derivable self-energies 



One expects that there should be a connection between the condition for particle conservation 
as expressed by A J = in Eq. (70), and the concept of a conserving approximation in the 



^-derivable sense. Below we show that AI = is always obeyed when the self-energy is 
^-derivable. 

We start by noting that Eq. (65) holds for any pair G(l, 2), E[G(1, 2)] provided E is of 
the ^-derivable form. In particular Eq. (65) does not assume that G, E[G] fulfill a Dyson 
equation. Therefore, given any orthonormal, but not necessarily complete basis, and 
writing (2(1,2) = . (^(r^Gj^Ti, r 2 )0*(r 2 ) we get from Eq. (65) after integrating over r 1; 

£ / dr'p^r, t')G^(t', t+) - G y (r- r')S^(r', r)] = 0, (74) 

which in matrix notation takes the form 

J dr'Tr[E(r, r')G(r', r+) - G(t~ , r')E(r', r)] = 0. (75) 

Here Ejj is the self-energy matrix obtained when the diagrams are evaluated using Gij and 
the from Eq. (3). The left hand side of Eq. (75), which is always zero for a ^-derivable 
E, can be written as Tr[A < (t, t)\ when A is given by Eq. (92) with B = E and C = G. It 
then follows from the general result (93) and the condition (73) that current conservation in 
the sense II = —Ir is always obeyed when E is ^-derivable. 

The above derivation of Eq. (75) relied on all the Coulomb matrix elements, V^i, being 
included in the evaluation of E. Thus the proof does not carry through if a general truncation 
scheme for the interaction matrix is used. However, in the special case of a truncated inter- 
action of the form (43), i.e. when the interaction is a two-point function, Eq. (75) remains 
valid. To show this, it is more appropriate to work entirely in the matrix representation and 
thus define $[Gij(r, r')] as the sum of a set of skeleton diagrams evaluated directly in terms 
Gij and V^-. With the same argument as used in Eq. (64), it follows that <3> is invariant under 
the transformation 

G^r, t') - e^G^r, r')e-^'\ (76) 

where A is now a discrete vector. By adapting the arguments following Eq. (64) to the discrete 
case we arrive at Eq. (65) with the replacements ri — > i and r 2 — > j and with the integral 
replaced by a discrete sum over j. Summing also over % leads directly to Eq. (75) which is the 
desired result. 

To summarize, we have shown that particle conservation in the sense II = —Ir, is obeyed 
whenever a ^-derivable self-energy is used and either (i) all Coulomb matrix elements V^i 
or (ii) the truncated two-point interaction of Eq. (43), are included in the calculation of E. 

4.4.1 Example: Transport through a single level 

As an illustrative example we consider transport through an Anderson impurity level con- 
nected to wideband leads, i.e. the retarded lead self-energies are taken to be frequency- 
independent equal to iY. The Hamiltonian of the central region reads 



He — Sc^c + Un^ny. 



(77) 
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Figure 5: Current-voltage characteristic for an Anderson impurity level with T = 0.65, e c = 
—4, U — 4. At finite bias, the non-selfconsistent GoW^otGnF] approximation yields different 
currents at the left and right interfaces (AI ^ 0) which means that charge conservation is 
violated. Moreover, GoWo predicts significant negative differential conductance for V > 0.8. 
In contrast, the ^-derivable HF and GW self-energies both conserve charge. 



where e c is the non-interacting energy and U is the correlation energy. 

In Figure 5 we show the current-voltage curve calculated from Eq. (68) using the self- 
consistent HF, selfconsistent GW, and non-selfconsistent (NSC) GW approximations. The 
parameters used are given in the figure caption. In the NSC calculation (referred to as G W ) 
we use the self-consistent Hartree-Fock GF to evaluate the GW self-energy. With this NSC 
self-energy we solve the Dyson's equation to obtain a NSC Green's function which is used to 
calculate the current. This "one-shot" approach is not conserving, and as a result the currents 
calculated in the left and right leads from Eq. (67) are not guaranteed to coincide. From the 
lower panel of figure 5 it can be seen that the G W self-energy does indeed violate charge 
conservation, and moreover leads to unphysical negative differential conductance for V > 0.8. 
On the other hand the selfconsistent HF and GW approximations both conserve charge. 

At this point we mention that the role of self-consistency in GW calculations has been 
much debated in the literature where it has been argued that G W , with G being the DFT 
Green's function, produces better band structures and band gaps than self consistent GW. 
The present example clearly demonstrates that, regardless of the performance of GW for band 
structure calculations, self-consistency is fundamental in nonequilibrium situations. 



5 TWO-LEVEL MODEL 



In this section we apply the general formalism presented in the preceding sections to a generic 
two- level model of a molecular junction. In this model the molecule is represented by its 
highest occupied (HOMO) and lowest unoccupied (LUMO) orbitals and the leads are repre- 
sented by one-dimensional tight-binding chains. With the aim of identifying universal trends 
we compare Hartree, Hartree-Fock, and GW calculations for the spectrum and IV character- 
istic. Not surprisingly the Hartree and HF results show large systematic differences due to 
the self-interaction errors in the Hartree potential. More interestingly, the dynamic correla- 
tions can have a large impact on both the spectrum and IV leading to significant deviations 
between the GW and HF results - in particular at finite bias. 

In section 5.1 we introduce the two-level model. In section 5.2 we study how dynamic 
screening effects influence the equilibrium position of the HOMO and LUMO levels both in 
the case of weak and strong coupling between molecule and the leads. In section 5.3 we 
consider the nonequilibrium transport properties of the model and explain the features of the 
IV curves in terms of the variation of the HOMO and LUMO positions as a function of the 
bias voltage. 
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Figure 6: The two-level model used to describe the HOMO and LUMO levels of a molecule 
couped to metallic leads, (a) The one-particle hopping matrix elements, (b) The electron- 
electron interactions. The interactions can be divided into intra-molecule interactions (Uq 
and Uhl) an d metal-molecule interactions (U ext ). 



5.1 Hamiltonian 



The model consists of two electronic levels coupled to two semi-infinite ID tight binding chains 
with nearest neighbor hopping, see Figure 6. The levels represent the HOMO and LUMO 
states of a molecule and the TB chains represent metallic leads. Electron-electron interactions 
on the molecule, and between the molecule and the first site of the chains are included. The 
Hamiltonian of the two-level model reads 

H — hi + h r + h mo i + h coup + U mo i + U ext . (78) 

Notice that we use a notation different from the canonical (L, C, i?)-notation introduced in 
section 2. This is because of the requirement that all interactions must be contained in the 
central region. Due to the interactions between the molecule and first sites of the leads, this 
implies that the central region should at least comprise the molecule and the first two sites 
of the leads. We enumerate the sites of the TB leads from — oo to —1 (left lead), and from 1 
to oo (right lead). Thus hi reads 

-l 

h i= Yl *(4x c m<r + 4 +la c ia ) (79) 

i=-oo <r=T,i 

with a similar expression for h r . The non- interacting part of the molecule's Hamiltonian reads 

hmol = ^ S £<* d aa d c*o ( 80 ) 
a=H,L <J=1,l 

and the coupling is given by 

Koup = t hyb (c\ a d aa + d\ a c la + cL la d aa + d ] aa C- 1(J ) (81) 

ol=H,Lo-=1,1 

For clarity we use c-operators for the lead sites and rf-operators for the HOMO/LUMO levels 
of the molecule. The interactions are given by 

Umoi = U (h m n Hl +h L ^n Ll ) + U HL n H n L (82) 
Uext = U ext (nin H + hih L + + n_in L ) (83) 

where e.g. % = d H *dHi + d H , dni is the number operator of the HOMO level, and fi\ is the 
number operator of the first site of the right lead. We set the Fermi level to zero corresponding 
to half-filled bands. In general, we write £l = £h + A , i.e. we use the non-interacting energy 
gap as a free parameter. The occupation of the molecule can then be controlled by adjusting 

5.2 Renormalization of molecular levels by dynamic screening 

In the low-bias regime, the transport properties are to a large extent determined by the 
position of the HOMO and LUMO levels relative to the Fermi energy of the metal electrodes. 
For this reason we shall first consider how the HOMO and LUMO position are affected by the 
interactions in the zero bias limit. The material presented in this section is part of ongoing 
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Figure 7: As a molecule approaches a metal surface its HOMO-LUMO gap is reduced 
by image-charge formation in the metal. If the molecule-metal bond is sufficiently strong 
(chemisorption) dynamic charge transfer between molecule and metal can give rise to addi- 
tional reduction of the gap. These renormalization effects requires a dynamic, i.e. frequency 
dependent, self-energy and thus cannot be described within the single-particle picture. 



work which will be published elsewhere. 

When a molecule is brought into contact with a metal surface a number of different 
mechanisms will affect the position of the molecule's energy levels. First, the levels are 
shifted by the electrostatic potential outside the surface. Second, hybridization effects shift 
and broaden the levels into resonances with finite life-times. For our model, the resonance 
width due to coupling to the leads is T pa \th y b\ 2 /t. Both the shift due to the surface potential 
and the hybridization are single-particle effects which can be described at a mean-field level 
such as Kohn-Sham (KS) or Hartree-Fock (HF) theory. On the other hand, correlation effects 
can renormalize the molecular spectrum in a way that cannot be described within a single- 
particle picture. Correlation effects are generally small for isolated molecules where HF usually 
yields good spectra, however, they can become significant when the molecule is in contact 
with a surface. An important example is the Kondo effect, where electronic interactions 
qualitatively changes the molecule's spectrum by introducing a narrow peak at the metal 
Fermi level (Goldhaber 1998 et al, Costi et al. 1994, Thygesen and Rubio et al. 2008). 
In weakly correlated systems such as molecules with a closed shell structure adsorbed at a 
surface, the effect of electronic interactions is expected to be less dramatic. However, as we 
shall see below, dynamic screening at the molecule-metal interface can introduce significant 
reductions of the HOMO-LUMO gap, which in turn will influence the transport through the 
molecule. Such screening effects can be observed in photoemission- and electron tunneling 
spectroscopy (Johnson and Hulbert 1987, Kubatkin et al. 2003, Repp et al. 2005). Recently, 
first-principles GW calculations for benzene physisorbed on graphite showed a HOMO-LUMO 
gap reduction of more than 3 eV due to substrate polarization (Neaton et al. 2006). More 
empirical treatments of polarization/screening effects using a scissors operator on the DFT 
spectrum have recently been applied to transport in molecules (Quek et. al 2007 and Mowbray 
et al. 2008) and scanning tunneling microscopy simulations (Dubois et al. 2007). So far the 
theoretical studies have focused on weakly coupled (physisorbed) molecules where the gap 



renormalization is induced by substrate polarization. This is the situation studied in section 
5.2.2 below. In section 5.2.3 we consider the to the case of strongly coupled (chemisorbed) 
molecules. 



5.2.1 The quasi-particle picture 

In interacting many-electron systems the concept of a single-particle eigenenergy becomes 
meaningless. However, for weakly correlated systems the concept can still be maintained due 
to the long life-time of certain states of the form cjj^o) (for 4> m unoccupied) or c m \^o) (f° r 
4> m occupied). These quasiparticle (QP) states describes the many-body iV-particle ground- 
state with and added electron (hole). The energy of the QP states relative to the groundstate 
energy, is given by the spectral function, A m (e) = — Im.G r rnm {e), where G r is the retarded 
Green's function. For weakly correlated systems A m (e) will be peaked at the QP energy, e m , 
which is equivalent to saying that the QP has a long life-time. It is instructive to notice that 
the QP energies measures the cost of removing/adding an electron to the state \4> m ) in the 
presence of interactions with the other electrons of the system. 

For non-interacting electrons the peaks in the spectral function coincide with the eigen- 
values of the single-particle Hamiltonian. Meanfield theories like KS or HF includes inter- 
actions at a static level, i.e. the single-particle eigenvalues correspond to the energy cost of 
adding/removing an electron when the state of all other electrons is kept fixed. This is the 
content of Koopman's theorem which states that (for <p m occupied), = E[^ F ]—E[c m ^ F ], 
i.e. the HF eigenvalues correspond to unrelaxed removal /addition energies. In general, the 
other electrons will respond to the added electron/hole and this will shift, or renormalize, 
the HF energies. The size of this effect is expected to qualitatively follow the linear response 
function, x, which gives the change in the particle density when the system is subject to an 
external field, 



This suggests a direct relation between the impact of dynamic relaxations, or screening, on 
the QP spectrum, and the static response function x(u — 0). 

5.2.2 Weak molecule-lead coupling 

We first consider the case of a weakly coupled, or physisorbed, molecule corresponding to 
small t hyb . We use the following default parameter values: t = 10, U = 4, U HL = 3, U ext = 2, 
thyb = 0.3, Ao = 4, which yield a resonance width of T ~ 0.01. In Fig. 8(left) we show the 
HOMO and LUMO positions as function of the interaction U ex t as calculated using the HF 
and GW approximations. As U ext is increased corresponding to the molecule approaching the 
surface, the GW gap decreases while the HF gap remains unchanged. In the simplest picture 
the gap reduction is due to the interaction between the added/removed electron and its image 
charge in the metal. This effect is not present in the HF single-particle spectrum: According 
to Koopman's theorem the added/removed electron interacts with the HF groundstate of the 
neutral system which contains no image charge. For small thyb where the molecule's levels 
are well defined in energy and localized it is possible to include the response of the metal 
to the added/removed electron by performing HF total energy calculation with constrained 




(84) 



HOMO/LUMO occupation numbers. Denoting by $q the minimizing Slater determinant 
with the constrain of n excess electrons on the molecule, we can define constrained HF energy 
levels as the total energy difference 

el% = HE\^ 1 ]-m])- (85) 

From the very good agreement between e^ f and e^ w seen in Fig. 8b, we conclude that GW 
includes the effect of relaxation or screening in the metal at the HF level. The situation is 
sketched in Fig. 10a. 
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Figure 8: Position of the HOMO and LUMO levels as function of interaction strength U ext 
(left) and non-interacting gap Ao (right) for a weakly coupled molecule (small T). The GW 
HOMO (LUMO) correspond to the HF energy cost of removing (adding) an electron to the 
molecule when the Slater determinant of the metal is allowed to relax. The gap reduction 
from (unrelaxed) HF to GW is thus due to polarization of the metal. This gap reduction is 
independent of A . 



Koopman's theorem allow us to write the difference between the HF single-particle lev- 
els and the result of constrained total energies, as (in case of the LUMO) ef F — = 
-Ef^o" 1 ] — -E[g^$q]. This energy difference has two contributions: A positive one from the 
interaction between the added electron and the induced density in the metal, and a negative 
one being the cost of forming the induced density. The classical image charge approximation 
in contrast assumes perfect screening in the metal and zero energy cost of polarizing the metal. 

From the right panel of Fig. 8 it can be seen that the renormalization of the gap is 
independent on the intrinsic gap of the molecule. This is expected since the image-charge and 
its interaction with the added electron/hole is independent on the HOMO-LUMO positions. 

According to the above, the size of the gap reduction for fixed U ext should depend on the 
polarizability of the metal. In Fig. 9(top) we show the dependence of the levels as function of 
t, i.e. the bandwidth of the metal. The effect is larger for small t corresponding to a narrow 
band. This is easily undestood by noting that narrow bands have larger DOS at Ef which in 
turn implies a larger density response function. Indeed, the right panel of Fig. 9 shows the 
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Figure 9: Left panels: Position of the molecule's HOMO and LUMO levels as a function of 
the metal bandwidth, t and the hybridization strength, thyb, respectively. GW (mol) refers to 
a calculation where only the interactions internally on the molecule have been treated within 
GW. Right panels: Static linear response functions (RPA) (<f) n \x{uJ = O)|0 n ) for the HOMO, 
LUMO and terminal site of the TB chain. The reduction of the correlated gap relative to HF 
is due to polarization of the metal and, for large th y b, of the molecule itself. Default parameter 
values are the same as in Fig. 8 



diagonal elements of the static (RPA) response function for the HOMO, LUMO and terminal 
site of the chain. The response function of the HOMO and LUMO is negligible for all values 
of t, while the response of the terminal site is significant and increases as t is reduced. 

5.2.3 Strong molecule-lead coupling 

We now turn to the case of a strongly coupled, or chemisorbed, molecule corresponding to 
non- negligible t hyb . In the bottom panel of Fig. 9 we show the center of the molecular 
resonances as a function of thyb- I n addition to the HF and GW values, we also show the 
result when only U mo i is treated at the GW level while U ex t is treated within HF. This allows 
us to isolate the correlation effects induced by the intra-molecular interactions from those of 
the metal-molecule interactions. Clearly, the correlated gap decrease relative to the HF gap as 
thyb is increased. It is also clear that the coupling-dependent part of the gap reduction comes 
from the interactions internally on the molecule, while the reduction due to V ext is largely 
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Figure 10: (al) Groundstate of a physisorbed molecule at a metal surface. (a2) When an 
electron is added to the molecule, the metal is polarized, (bl) Groundstate of a chemisorbed 
molecule at a metal surface. (b2) When an electron is added to the LUMO the metal is 
polarized, and charge is transfered from the molecule to the metal. The screening effects 
stabilize the system with added/removed electron and this shift the occupied (unoccupied) 
quasi-particle levels up (down). 



independent of thyb- Since V mo i does not produce any renormalization of the levels of the 
free molecule (see the t^yb — ► limit), the mechanism responsible for the gap reduction must 
involve the metal. From the lower right panel of Fig. 9, we see that the response functions 
of the HOMO and LUMO states increases with thyb indicating the gap reduction due to V mo i 
is of a similar nature as the image charge effect, but with the molecule itself being polarized. 
The effect increases with thyb because charge transfer between the molecule and the metal due 
to the external field from the added/removed electron, is larger when resonances are broad 
and have larger overlap with the metal Fermi level. The situation is sketched in Fig. 10b. 

5.3 Nonequilibrium transport 

The analysis of the previous sections show that dynamic screening effects can have a large 
effect on the spectrum of the molecule in contact with leads. In this section we shall see 
that the application of a bias voltage leads to additional renormalization of the spectrum. 
For simplicity we limit the model to include intra-molecule interactions, i.e. we set U ex t to 
zero. This means that reduction of the HOMO-LUMO gap due to image charge formation in 
the leads is not included. Whereas the presence of intra-molecule interactions did not have a 
large influence on the equilibrium positions of the HOMO and LUMO levels for small values 
of T (see the GW (mol) result in Fig. 9), we will see that this is no longer true under finite 
bias conditions, where intra-molecular screening is strongly enhanced and the life-times of the 
HOMO and LUMO levels can be significantly reduced due to QP scattering. Both effects lead 



to a reduction of the HOMO-LUMO gap as function of the bias voltage with a large impact 
in the calculated IV curve. 
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Figure 11: The HOMO and LUMO resonances of the two-level model under zero and non-zero 
bias voltage. As indicated the bias affects both the position and width of the resonances and 
this will in turn affect the dl/dV curve. Reprinted with permission from Thygesen, Phys. 
Rev. Lett. 100, 166804 (2008). Copyright 2008 by the American Physical Society. 

Throughout this section we use following parameters: A = 2, U — 2, Uhl = 1-5, 
t = 10. By varying the one-particle energy £ , we can control the equilibrium occupation 
of the molecule, N e \. We consider the case of weak charge transfer to the molecule, i.e. N e \ 
ranges from 2.0 to 2.1, corresponding to Bp lying in the middle of the gap and slightly below 
the LUMO, respectively. The Fermi level is set to zero, and the bias is applied symmetrically, 
i.e. /jLl = V/2 and /j,r = —V/2. The situation is illustrated in figure 11. 

In Fig. 12 we show the calculated dl/dV curves (obtained by numerical differentiation of 
I(V)) for different values of V and N e \. We first notice that the 2B and GW approximations 
yield similar results in all the cases indicating that the higher order terms in the GW self- 
energy are fairly small. For T = 1.0, all methods yield qualitatively the same result. For 
even larger values of T (not shown), and independently of N e \, the results become even more 
similar. In this strong coupling limit, single-particle hybridization effects will dominate over 
the interactions and xc-effects are small. 

Focusing on the T = 0.25 case we see that the Hartree approximation severely overesti- 
mates the low-bias conductance. This is a consequence of the self-interactions (SI) contained 
in the Hartree potential which leads to an underestimation of the (equilibrium) HOMO-LUMO 
gap, see figure 13 for V = 0. On the other hand the HF, 2B, and GW methods lead to very 
similar conductances in the low-bias regime. This is consistent with the results of the previous 
section which showed that intra-molecular correlations do not renormalize the equilibrium HF 
gap much for small T. For T = 1.0 the slightly larger conductance in GW and 2B is due to 
the slight reduction of the equilibrium gap. 

We notice that the lower left graph (r = 0.25, iVei = 2.0) shows an interesting feature. 




Figure 12: dl/dV curves for different values of the tunneling strength T and occupation of the 
molecule, N e \. The curves are calculated using different approximations for the xc self-energy. 
Reprinted with permission from Thygesen, Phys. Rev. Lett. 100, 166804 (2008). Copyright 
2008 by the American Physical Society. 



Namely, the HF, 2B, and GW curves all contain an anomalously strong conductance peak. 
Interestingly, the peak height is significantly larger than 1 which is the maximum conductance 
for a single level (the Anderson impurity model). Moreover, the full width at half maximum 
(FWHM) of the peak is only ohf = 0.27 and ct 2 b/gw — 0.12, respectively, which is much 
smaller than the tunneling broadening of 2r = 0.5. We note in passing that the peak looses 
intensity as N e \ is increased, and that the Hartree approximation does not reproduce the 
narrow peak. 

5.3.1 Influence of bias on the HOMO and LUMO positions 

To understand the origin of the anomalous conductance peaks, we consider the evolution of the 
HOMO and LUMO positions as a function of the bias voltage, see Fig. 13 (the 2B result is left 
out as it is similar to GW). Focusing on the upper panel of the figure (corresponding to N d = 
2.0), we notice a qualitative difference between the Hartree and the Si-free approximations: 
While the Hartree gap expands as the levels move into the bias window, the HF and GW gaps 
shrink leading to a dramatic increase in current around V = 2.5 and V = 1.3, respectively. 
This is clearly the origin of the anomalous dl/dV peaks. But why do the Si-free gaps shrink 
as the bias is raised? 
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Figure 13: Position of the HOMO and LUMO levels as a function of the bias voltage for 
the Hartree (crosses), HF (triangles), and GW (circles) approximations. The horizontal lines 
show the FWHM of the GW resonances. The FWHM of the Hartree and HF resonances is 2r 
independently of V. Notice the differences in the way the levels enter the bias window: The 
Hartree gap opens while the HF and GW gaps close. In the upper graph T = 0.25, N c \ = 2.0 
(symmetric case). In the lower graph T = 0.25, N e \ = 2.1. Reprinted with permission from 
Thygesen, Phys. Rev. Lett. 100, 166804 (2008). Copyright 2008 by the American Physical 
Society. 



Let us consider the change in the HOMO and LUMO positions when the bias V is in- 
creased by 25V. In general this change must be determined self-consistently, however, a 
"first iteration" estimate yields a change in the HOMO and LUMO occupations of 8nn ~ 
— Ah(— V/2)SV and 5til ~ Al(V/2)8V, respectively. Here A H / L is the spectral function, or 



equivalently the DOS, of the HOMO/LUMO levels. At the HF level this leads to 



Se H « [-U A(-V/2) + 2U HL A(V/2)]5V (86) 
5e L « [C/q^/2) - 2C/ J/L A(-y/2)]5y (87) 

where A = Hh + is the total DOS on the molecule and we have used that Ah(—V/2) 
A(— V/2) and A(V/2) w A(V/2). The factor 2 in front of Uhl accounts for interactions with 
both spin channels. In the symmetric case (N A = 2.0) we have A(—V/2) = A(V/2). Since 
U < 2U HL this means that Se H > and 5e L < 0, i.e. the gap is reduced as V is raised. 
Moreover it follows that the gap reduction is largest when A(±V/2) is largest, that is, just 
when the levels cross the bias window. In the general case (N e \ ^ 2.0) the direction of the 
shift depends on the relative magnitude of the DOS at the two bias window edges: a level 
will follow the edge of the bias window if the other level does not intersect the bias edge. It 
will move opposite to the bias, i.e. into the bias window, if the other level is close to the bias 
window edge. This is effect clearly seen in the lower graph of Fig. 13 (triangles). Thus the gap 
closing mechanism has the largest impact on the dl/dV curve when the HOMO and LUMO 
levels hit the bias window simultaneously. Moreover, the effect is stronger the larger Uhl/Uq, 
and the smaller T (the maximum in the DOS is ~ 1/T). At the Hartree level, Eqs. (86) and 
(87) are modified by replacing U by 2U due to self-interaction. This leads to an effective 
pinning of the levels to the bias window which tends to open the gap as V is increased, see 
Fig. 13 (crosses). 

The above analysis shows why the HF gap is reduced as the levels hit the bias window. 
Interestingly, the bias-driven gap reduction is even stronger in GW and as a consequence the 
GW conductance peak occurs at much lower bias (V = 1.5) than the HF peak (V = 2.5). 
Part of the downshift of the GW conductance peak can be explained from the smaller GW 
equilibrium gap. Indeed, for V = the HF gap is ~ 0.3 larger than the GW gap. However, 
this effect alone cannot account for the large down-shift. 

The reason for the bias induced reduction of the GW gap is two-fold: First, intra-molecular 
screening effects are enhanced as the chemical potentials move closer to the molecular levels 
and increase the susceptibility of the levels. This is analogue to the (equilibrium) situation 
of increasing T shown in lower panels of figure 9. The susceptibility of a molecular level is 
roughly given by the magnitude of the level's DOS at E F (or chemical potentials). In the 
latter case this is achieved by broadening the resonances; in the former case by bringing the 
chemical potential(s) closer to the levels. Secondly, the rate of QP scattering, i.e. the rate 
at which the initial state cj |^o) is destroyed due to electron-electron interactions, increases 
with the bias. This follows from Fermi's Golden rule by realizing that the number of available 
final states of the form Cj|^o) having the same energy as the initial state, increases with bias. 
The enhanced QP scattering reduces the life-time of the HOMO and LUMO QPs, which is 
equivalent to a broadening of the molecular resonances. 

The full width at half maximum (FWHM) of the GW resonances is indicated by horizontal 
lines in Fig. 13. For low bias, the width of the GW resonances is the same as the width of 
the Hartree and HF resonances. The latter is determined by the coupling to leads and equals 
2r = 0.5 independently of the bias. According to Fermi-liquid theory, QP scattering at the 
Fermi level is strongly suppressed in the ground state, i.e. ImEj^ep) = for V = (recall 



that ImS is inversely proportinal to the life time). However, as the bias is raised the phase 
space available for QP scattering is enlarged and ImE increases accordingly. As a result of the 
additional level broadening, A(±V/2) increases more rapidly as a function of V. According 
to Eqs. (86,87), this will accelerate the gap closing reduction already at the HF level. Finally, 
we notice that the long, flat tails seen in the dl/dV of the GW /2B calculations are also a 
result of the spectral broadening due to QP scattering. 



In this section we combine the GW scheme with a Wannier function (WF) basis set to study 
electron transport through two prototypical junctions, namely a benzene molecule coupled to 
featureless leads and a hydrogen molecule between two semi-infinite Pt chains. In section 6.1 
we briefly present the computational scheme. In the following two sections we analyze the 
energy spectrum and transport properties of the benzene junction. Finally in section 6.4 we 
present results for the IV curve of the Pt-H 2 -Pt junction. 

6.1 Computational details 

Below we review our computational scheme for GW transport calculations in a WF basis 
discussed in detail in Thygesen and Rubio 2008. In a first step, periodic supercell DFT calcu- 
lations are performed for the leads as well as the central region containing the molecule plus 
part of the leads. We use the Dacapo code (Bahn and Jacobsen 2002) which applies ultrasoft 
pseudopotentials (Vanderbilt 1990) for the ion cores. The KS eigenstates are expanded in 
plane waves with a cut off energy of 340 Ry and the PBE xc-functional (Perdew et al. 1996) 
is used. In the second step, the KS eigenstates are transformed into maximally localized, 
partially occupied WFs, {0 n (r)}, and the KS Hamiltonians of the central region and the 
leads are subsequently evaluated in terms of the WF basis. The eighteen maximally localized 
WFs obtained for the benzene molecule are shown in Fig. 14. The xc-potential is excluded 
from the Hamiltonian of the central region in order to avoid double counting when the GW 
self-energy is added. The central region Hamiltonian reads 



where v ps is the pseudopotential and Vh is the Hartree potential. Notice that v ps and Vh con- 
tain contributions from the ion cores and electron density of the leads as they are obtained 
from a supercell calculation with part of the leads included. 

The Coulomb integrals are evaluated for the WFs of the central region, 



For the correlation part of the GW self-energy, S corr = T, GW — E x , we use the effective inter- 
action introduced in section 3.1, i.e. only Coulomb integrals of the form U^ and Vijji are 
included. For the Hartree and exchange self-energies, and S x , which are easily evaluated 
from Eqs. (58,59), we use all the Coulomb matrix elements. Notice, that we need even 
though the Hartree potential from electrons in C is already contained in tv The reason is 



6 PROTOTYPE MOLECULAR JUNCTIONS: C 6 H 6 and H 2 



[h c ]ij = (<fii\ - ^ v2 + v vs + Vh\<frj), 



(88) 




(89) 



that the latter is the equilibrium Hartree potential of the DFT calculation, which might well 
differ from the Hartree potential of a nonequilibrium GW calculation. 

The retarded Green's function is evaluated from 

G r = [u _ hc _ Ei _ Er _ (E r [G ] _ EJ^M]) - Kc[G\\-\ (90) 

where the frequency dependence has been omitted for notational simplicity. Several comments 
are in order. First, £l and are the lead self-energies of Eq. (41) (in the wide-band 
approximation and £r are diagonal and frequency independent). The term Avh = T> r h [G} — 
Yj^lgi^] is the change in Hartree potential relative to the equilibrium DFT value. This change 
is due to the applied bias and the replacement of v xc by Y> xc . In this work Y> xc can be either 
the exchange or the GW self-energy. Finally we notice that the bias dependence of the various 
quantities entering Eq. (90) has been suppressed for notational simplicity. 

a) b) 




Figure 14: (a) Illustration of a benzene molecule coupled to featureless electrodes with dif- 
ferent chemical potentials, (b) Iso-surfaces for the 18 partially occupied Wannier functions 
used as basis functions in the calculations. The WFs are linear combinations of Kohn-Sham 
eigenstates obtained from a DFT-PBE plane-wave calculation. 



6.2 Equilibrium spectrum of benzene 

In Fig. 15 we show the total density of states (DOS) of the isolated benzene molecule calculated 
using three different approximations: (i) DFT-PBE (ii) Hartree-Fock (iii) fully self-consistent 
GW. The DOS is given by 

D(e) = —J2^GUe), (91) 

n=l 

where the sum runs over all WFs on the molecule, and the GF is obtained from Eq. (90) 
using a wide-band lead self-energy of T = 0.05. We stress that our calculations include the 
full dynamic dependence of the GW self-energy as well as all off-diagonal elements. Thus no 
analytical extension is performed, and we do not linearize the self-energy around the DFT 
eigenvalues to obtain an approximate quasi-particle equation as is done in standard GW cal- 
culations. 
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Figure 15: Density of states for a benzene molecule weakly coupled to featureless leads (T = 
0.05). The common Fermi levels of the leads is indicated. Notice the characteristic opening of 
the band gap when going from DFT-PBE to HF, and the subsequent (slight) reduction when 
correlations are included at the GW level. Reprinted with permission from Thygesen and 
Rubio, Phys. Rev. B 77, 115333 (2008). Copyright 2008 by the American Physical Society. 



The spectral peaks seen in Fig. 15 occurring above (below) the Fermi level correspond to 
electron addition (removal) energies. In particular, the HOMO level should coincide with the 
ionization potential of the isolated molecule, which in the case of benzene is J cxp = —9.2 eV 
(NIST Chemistry WebBook). The PBE functional overestimates this value by 3 eV, giving 
/pbe = —6.2 eV in good agreement with previous calculations (Niehaus et al. 2005). The 
HF and GW calculation yields Irf = —9.7 eV and I GW = —9.3 eV, respectively. Given the 
limited size of the Wannier basis, the precise values should not be taken too strict. However, 
the results demonstrate the general trend: KS theory with a local xc-functional underestimates 
the HOMO-LUMO gap significantly due to SI errors; HF overestimates the gap slightly; GW 
reduces the HF gap slightly through the inclusion of dynamic screening. 

In Fig. 16 we plot the size of the HOMO-LUMO gap as a function of the coupling strength 
T. The position of the levels has been defined as the first maximum in the DOS to the left 
and right of the Fermi level. Both the HF and GW gaps decrease as T is increased. This 
observation is consistent with the model calculations of section 5.2.3 where it was found that 
the gap reduction due to the correlation part of the GW self-energy, A corr , can be understood 
as a virtual charge transfer between molecule and leads. The reduction of the HF gap as 
function of T is a consequence of the redistribution of charge from the HOMO to the LUMO 
when the resonances broaden and their tails start to cross the Fermi level. This is completely 
analogue to the bias induced gap reduction discussed in section 5.3. These results for benzene 
show that the conclusions obtained from the two-level model apply to more realistic systems. 
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Figure 16: The HF and GW HOMO-LUMO gap of the benzene molecule as a function of 
the coupling strength V. The difference between the curves represents the reduction in the 
gap due to the correlation part of the GW self-energy. This value increases with the coupling 
strength, as polarization of the molecule via dynamic charge transfer to the metal becomes 
possible (see section 5.2.3). Reprinted with permission from Thygesen and Rubio, Phys. Rev. 
B 77, 115333 (2008). Copyright 2008 by the American Physical Society. 



6.3 Conductance of benzene 

We consider electron transport through the benzene junction under a symmetric bias, [il/r = 
±V/2, and a wide-band coupling strength of Tl = Tr = 0.25 eV. In Fig. 17 we compare the 
differential conductance, dl/dV, calculated from self-consistent DFT-PBE, HF, and GW, as 
well as non self-consistent GqWq using either the DFT-PBE or HF Green's function as Go- 
The dl/dV has been obtained by numerical differentiation of the I(V) curves calculated from 
Eq. (68). For the DFT calculation the finite-bias effects have been included at the Hartree 
level, i.e. changes in the xc-potential have been neglected. We notice that the HF and 
GoHo^hf] results are close to the self-consistent GW result. These approximations all yield 
a nearly linear IV with a conductance of ~ 0.05Go- in contrast the DFT and GoHo[Gt>ft] 
yield significantly larger conductances which increase with the bias voltage. We note that the 
violation of charge conservation in the GqWq calculations is not too large in the present case 
(A/// < 5%). This is in line with our general observation, e.g. from the Anderson model, 
that AI/I increases with I. 

The trends in conductance can be understood by considering the (equilibrium) DOS of the 
junction shown in Fig. 18. As for the weakly-coupled (free) benzene molecule whose spectrum 
is shown in Fig. 15, the DFT HOMO-LUMO gap is much smaller than the HF gap, and this 
explains the larger DFT conductance. The GW gap falls in between the DFT and HF gaps, 
however, the magnitude of the DOS at Ep is very similar in GW and HF which is the reason 
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Figure 17: Differential conductance of the benzene junction for = = 0.25 eV. Notice 
the strong Go dependence of the GqWo result. Reprinted with permission from Thygesen and 
Rubio, Phys. Rev. B 77, 115333 (2008). Copyright 2008 by the American Physical Society. 



for the similar conductances. It is interesting to notice that the HOMO-LUMO gap obtained 
in the G W calculations resemble the gap obtained from G , and that the self-consistent GW 
gap lies in between the Go^o[Gt>ft] and Go^o^hf] gaps. 

The increase in the GoVFoIGdft] conductance as a function of bias occurs because the 
LUMO of the G W / o[G'dft] calculation moves downwards into the bias window and becomes 
partly filled as the voltage is raised. In a self-consistent calculation this would lead to an 
increase in Hartree potential which would in turn raise the energy of the level. The latter 
effect is missing in the perturbative GqWq approach and this can lead to uncontrolled changes 
in the occupations as the present example shows. 

Finally we notice that the GoWo[Ctdft] DOS is significantly more broadened than both 
the GoWo[Chf] and GW DOS. The reason for this is that, as a direct consequence of the 
small HOMO-LUMO gap, DFT yields a the larger DOS close to E F . The larger DOS in 
turn enhances the QP scattering and leads to shorter life-times of the QP in the GoWo[Gt>ft] 
calculation. Since the QP life-time is inversely proportional to ImT.cw this explains the 
broadening of the spectrum. 

6.4 Pt-H 2 -Pt junction 

We consider a molecular hydrogen bridge between infinite atomic Pt chains as shown in the 
inset of Fig. 19. Experimentally, the conductance of the hydrogen junction is found to be 




Figure 18: Equilibrium DOS for the benzene molecule coupled to wide-band leads with a 
coupling strength of Fl — = 0.25 eV. Upper panel shows DFT-PBE and HF single- 
particle approximations while the lower panel shows the self-consistent GW result as well as 
one-shot GoWo results based on the DFT and HF Green's functions, respectively. Reprinted 
with permission from Thygesen and Rubio, Phys. Rev. B 77, 115333 (2008). Copyright 2008 
by the American Physical Society. 



close to the conductance quantum, Go = 2e 2 /h (Smit et al. 2002), and this value has been 
reproduced by DFT calculations (Thygesen and Jacobsen 2005). Below we present CPU- 
transport results for a simplified model of this system (using infinite Pt chains as leads), and 
refer to Thygesen and Rubio 2007 for further details on the calculations. 

In the upper panel of Fig. 19 we show the local density of states (LDOS) at one of the two 
H orbitals as calculated within DFT using the PBE xc-functional, as well as self-consistent 
HF (in the central region). In DFT the H 2 bonding state is a bound state at —7.0 eV relative 
to E F , while the anti-bonding state lies at 0.4 eV and is strongly broadened by coupling to 
the Pt. Moving from DFT to HF the bonding state is shifted down by ~ 8 eV because for 
occupied states the exchange potential is more negative than the DFT xc-potential. The 
same effect tends to drive the half-filled anti-bonding state down but in this case the resulting 
increase in the Hartree potential (about 4 eV) stops it just below E F . 

In the lower panel of Fig. 19 we show the LDOS calculated in GW as well as GqWq starting 
from either DFT or HF, i.e. Go is either Gdft or Ghf- The large deviation between the two 
GqWq results is not surprising given the large difference between Gdft and Ghf- A closer 
analysis of the origin of this deviation can be found elsewhere (Thygesen and Rubio 2007). 




Figure 19: Local density of states at one of the H orbitals of the Pt-H-H-Pt contact shown 
in the inset. Reprinted with permission from K. S. Thygesen and A. Rubio, J. Chem. Phys. 
126, 091101 (2007). Copyright 2007, American Institute of Physics. 



1 

0.8 



0.6 

CVJ 






V(volt) 



Figure 20: I-V and dl/dV for the hydrogen contact as calculated in DFT(PBE) and self- 
consistent GW . V is the source-drain bias voltage. Reprinted with permission from K. S. 
Thygesen and A. Rubio, J. Chem. Phys. 126, 091101 (2007). Copyright 2007, American 
Institute of Physics. 



We are, however, aware that part of this large difference could be due to the limited size of 
the basis. We also mention that the LDOS results of Fig. 19 can be largely reproduced by 



including only the second-order GW diagram in the self-energy. Thus the higher-order RPA 
diagrams are less important in this case. 

In Fig. 20 we show the self-consistently calculated IV characteristics in DFT and GW. 
At low bias both schemes yield a conductance close to the experimental value of 1G - The 
DFT conductance is nearly constant over the bias range, and is in fact very similar to the 
HF result (not shown). In contrast the GW conductance falls off at larger bias. This is due 
to enhancement of quasi-particle scattering at finite bias. The QP scattering reduces the life 
time of the QPs leading to broadening of the spectral peak associated with the anti-bonding 
state of the hydrogen molecule in agreement with the results for he two level model discussed 
in section 5.3. Since the correlation induced life-time of QP at the Fermi level, ImS corr (£'j?), 
vanishes identically in equilibrium, the finite-bias conductance suppression seen in figure 20 
is a direct result of the non-equilibrium treatment of correlations. 

7 SUMMARY AND PERSPECTIVES 

The feasibility of using many-body perturbation theory in combination with the Keldysh 
Green's function formalism to address nonequilibrium quantum transport through an inter- 
acting region coupled to non- interacting leads has been demonstrated. The effect of elec- 
tronic correlations was incorporated into the Green's function via the GW self-energy, and 
the coupling to leads was treated exactly (to all orders in the hopping between leads and 
central region). The important connection between self-consistency in the GW self-energy 
and charge conservation was emphasized, and it was demonstrated that a non self-consistent 
treatment of the GW self-energy (the G W approach) violates the continuity equation and 
produce unphysical results at finite bias. This, together with the arbitrariness of the G W 
approximation due to its Go-dependence, speaks in favor of the self-consistent GW approach 
to nonequilibrium transport. 

The role of dynamic correlation effects in quantum transport was illustrated by applying 
the GW-transport scheme to a generic two-level model, a benzene molecule between feature- 
less leads, and a hydrogen molecule between infinite Pt chains. It was shown that dynamic 
polarization of the leads as well as the molecule itself can lead to significant reduction of the 
molecule's HOMO-LUMO gap. The polarization effects were found to increase with the bias 
voltage where also quasiparticle scattering is strongly enhanced leading to broadening of the 
molecular resonances. As was shown, all these effects can have a large impact on the calcu- 
lated IV curve. This should always be kept in mind, when interpreting results of meanfield 
(DFT or Hartree-Fock) transport calculations - in particular under finite bias. 

As mentioned in the introduction, the quantitative theoretical description of quantum 
transport in nano-scale structures from first principles is an extremely complex problem. 
Nevertheless, simulations methods with predictive power are required to advance the field 
further. It has been known for several years that the standard DFT-NEGF scheme fails to 
predict even the zero-bias conductance of certain classes of systems. This state of affairs 
makes it difficult, although not impossible, to link theory and experiments and thereby stim- 
ulate the development of nano-scale electronics. 



As illustrated by the examples given in this chapter, reliable schemes for quantum trans- 
port should account for dynamic correlation effects in some way or another. The GW method 
discussed here includes some correlation effects, but misses others, e.g. the side peaks in the 
spectral function of the Anderson model are not well reproduced (Thygesen and Rubio 2007). 
Methods developed for strongly correlated systems, such as density matrix renormalization 
group theory (Schmitteckert and Evers 2008), are limited to simple models due to their in- 
expedient scaling with system size. The effective single-particle scheme of TDDFT makes 
it an attractive alternative to many-body perturbation theory, in particular for dynamical 
transport phenomena (Stefanucci et al. 2007). However, the inclusion of correlation effects 
requires use of xc-potentials with memory which have so far proved difficult to construct. 

Another important aspect of the problem is related to the coupling between electrons and 
nuclei. Despite the large difference in the general time scales of electronic and nuclear motions, 
electronic wavepackets quite often couple with the dynamics of nuclear motion (Frederiksen et 
al. 2004, Verdozzi et al. 2006). The proper incorporation of the electronic-nuclear interaction 
is crucial for describing a host of dynamical processes including Joule heating, electromi- 
gration, laser-induced electronic transport and electron transfer in molecular, biological, or 
electrochemical systems. Within the groundstate DFT framework, the computation of forces 
on the nuclei is trivial thanks to the Hellman-Feynman theorem. The situation is more com- 
plex out of equilibrium, and even more so in combination with a many-body description of 
the electrons, where the Hellman-Feynman theorem does not apply. However, the electron- 
ion dynamics must eventually be taken properly into account for a realistic description of a 
large class of molecular devices relevant for technological applications such as fast, integrated, 
optoelectronic nanodevices. 
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APPENDIX 

Let B(t, t') and C{t, t') be two matrix valued functions on the Keldysh contour, and consider 
the commutator A defined by 



where matrix multiplication is implied. Under steady state conditions where the real time 
components of B and C can be assumed to depend only on the time difference t' — t, the 



Arina) . 




(92) 



following identity holds: 

Tr[A<(t,t)] = J ^Tr[B<(u)C>(u) - B>(u)C<(u)]. (93) 

To prove this relation we first use the Langreth rules to obtain 

A<(t,t') = J [S < (i,i 1 )C°(*i,0 + 5 r (t,ti)C < (ti,0 

- C<(t, *i)S a (t!, - C r (t, fi)B<(f i, 0] dti. 

Since all quantities on the right hand side depend only on the time difference we identify the 
integrals as convolutions which in turn become products when Fourier transformed. We thus 
have 



A<(t,t) = J^A<(u) 



I 



^-[B<(u)C a (u) + B r (u)C<(u) 

Z7T 



-C < (w)S°(w) - C r {uj)B < {uj)]. 
Eq. (93) now follows from the cyclic property of the trace and the identity G r — G a = G > — G < . 
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